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ABSTRACT 

Wc perform an autocorrelation study of the Auger data with the aim to constrain the number 
density n s of ultrahigh energy cosmic ray (UHECR) sources, estimating at the same time the effect 
on n s of the systematic energy scale uncertainty and of the distribution of UHECR. The use of global 
analysis has the advantage that no biases are introduced, cither in n s or in the related error bar, by 
the a priori choice of a single angular scale. The case of continuous, uniformly distributed sources is 
nominally disfavored at 99% C.L. and the fit improves if the sources follow the large-scale structure 
of matter in the universe. The best fit values obtained for the number density of proton sources are 
within a factor ~2 around n s ~ 1 x 10~ 4 /Mpc 3 and depend mainly on the Auger energy calibration 
scale, with lower densities being preferred if the current scale is correct. The data show no significant 
small-scale clustering on scales smaller than a few degrees. This might be interpreted as a signature of 
magnetic smearing of comparable size, comparable with the indication of a « 3° magnetic deflection 
coming from cross-correlation results. The effects on the above results of some approximations done 
is also discussed. 

Subject headings: cosmic rays — large-scale structure of universe — methods: statistical 



1. INTRODUCTION 

Evidence is now emerging that ultrahigh energy 
cosmic rays (UHECRs) have an astrophysical ori- 
gin, as opposed to being generated in exotic top- 
down models: The detection of a spectral suppres- 
sion consistent with the Greisen-Zatsepin-Kuzmin 
(GZK) effect (jGreisen '66l ; IZatsepin fc Kuzmin '6rf ) 
by b oth HiRes (lAbbasi et al. '08( 1 and 
Auger ((Abraham et al. '08d( l collaborations, to- 
gether with the Auger bounds o n the UHE 
neutrino flux (lAbraham et al. '08bD. on the 
photon fraction ((Abraham et al '08a( l and on 
the anisotropy tow ar ds the Galactic Cen- 
ter (|Aglietta et al '07l ; lAloisio fc Tortorici '07( l are 
all consistent with this scenario. The next step is 
clearly the identification of the sources of UHECRs, 
an arena where anisotropy studies play a crucial role. 
Yet, the limited angular resolution of extensive air 
shower detectors and especially the deflections that 
charged particles suffer in astrophysical magnetic fields 
make the task highly non trivial. This is especially 
troublesome given that the UHECRs chemical compo- 
sition is unknown, that we lack a detailed knowledge 
of the Galactic magnetic field structure and, above all, 
of the very magnitude and structure of extragalactic 
magnetic fields (EGMF) outside of cluster cores. These 
limitations — together with the small statistics available 
at present — suggest that, at least in an initial phase, 
charged particle astronomy may be limited to the in- 
ference on the statistical properties of UHECR sources, 
rather tha n a detailed study of single accelerators. 

In Ref. ((Cuoco et al. '08bh . we found that a global 
comparison of the two-point auto-correlation function 



of the data with the one of catalogues of potential 
sources is a powerful diagnostic tool: This observable is 
less sensitive to unknown deflections in magnetic fields 
than the cross-correlation function, while keeping a 
strong discriminating power among source candidates. 
In particular, the autocorrelation function of (sub-) 
classes of galaxies have different biases with respect to 
the large-scale structure (LSS) of matter. As a result, 
the best fit value for the density n s of different source 
classes may differ, especially if only one or a small 
range of angular scales is considered. Although the bias 
of different source classes differs maximally at small 
angular scales, we showed that the statistically most 
significant differences are at intermediate angular scales, 
where both the larger number of cosmic ray pairs (CR) 
and of galaxy pairs leads to relatively smaller error bars. 
Moreover, the autocorrelation function on larger angular 
scales becomes less dependent on possible deflections in 
the Galactic and extragalactic fields. 

In this article we derive the number density of 
UHECR sources using the recently published ar- 
rival di rections and energ ies of the 27 Auger 
events ((Abraham et al. '08c( l with estimated en- 
ergy E > 57EeV, thereby complementing the 
study ((Cuoco et al. '08a ) with a concrete example 
for a comparison of the global cumulative autocorre- 
lation function of sources and UHE CRs. Note that 
we showed in Ref. ((Cuoco et al. '08bT ) that, even in an 
idealized case where systematics play no major role, 
roughly three times the numb er of "useful" events that 
can be extracted from Ref. (|Abraham et al. '08d ) are 
required to start distinguishing between different sub- 
classes of sources. Thus a study of the kind envisaged 
in Ref. ((C uoco et al. '08bf l is unrealistic at present. 
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Here, we restrict ourselves to more modest goals: i) To 
compare the data against predictions of two toy model 
cases of uniformly distributed sources and of "normal 
galaxies" (i.e. source s that have the sam e distribution as 
the PSCz catalogue ([Saunders et al. '"00l )) which we shall 
refer to with the two values for the label k = {uni, LSS}, 
respectively, ii) To study the effect on the allowed range 
of n s of a systematic error on the energy scale of the 
UHECR experiment. Note that preliminary results of 
the clustering o f the Auger events has been presented in 
(|Mollerach '07ft , but astrophysical implications have not 
been discussed there. 

We review the statistical method we use in Sec. [2j and 
apply it in Sec. [3] to the Auger data, providing some in- 
terpretation of the results. In Sec. [4] we discuss some 
limitations and caveats of the analysis. Finally, wc sum- 
marize in Sec. [5j 

2. STATISTICAL ANALYSIS OF THE DATA 

The use of correlation functions is well suited to the 
study of over- and underdensities of non-uniformly dis- 
tributed sources and of the resulting anisotropies in the 
radiation received from them. Since the number of CR 
events published by Au ger is still small, we use in our 
analysis following Ref. ( Kachelriefi fc Semikoz '06ft the 
cumulative two-point autocorrelation function C($) de- 
fined as 
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i.e. the number of pairs within the angular distance t9. 
Here, N is the number of CRs considered, $y is the 
angular distance between events i and j and O is the 
step function (with 0(0) = 1). 

This function is straightforward to compute for the 
data, and denoted then by C* ($). For a specific model 
hypothesis X, a set of functions d(fl\X) is obtained in 
the following way: Sources with equal luminosities are 
distributed within a sphere of 180/i _1 Mpc either uni- 
formly or following t he three-dimensional LSS as given by 
a smoothed version (ICuoco et al. '061: ICuoco et al. '08ah 
of the PSCz catalogue ([Saunders et al. '00ft . For the 
LSS case (but not for the uniform case) sources and CR 
events within the PSCz mask are excluded, leaving 22 
CR events. Note that the mask mostly overlaps with the 
Galactic plane and bulge region, where larger deflections 
due to the Galactic magnetic field are expected: The 
mask is thus not only a catalogue limitation, but also 
implements to some extent a physically motivated angu- 
lar cut. Finally, each source k, at redshift Zk, is weighted 
by the factor 

1 f°° 

-j / E- s dE (2) 

z k J Ei(E oxlt ,z k ) 

accounting for its redshift dependent flux suppression 
and the CR energy losses. These are parametrized as 
a continuous process in term of the function Ei(Ef,Zk) 
which gives the initial injection energy Ei for a parti- 
cle leaving the source at Zk and arriving at the Earth 
wit h energy Ef < Ei. Further details are given 
in (|Cuoco et al. '06ft . The injection spectral index s is 
assumed to be the same for all the sources and equal 
to 2.0. The dependence on s is however weak as shown 



in more detail in (|Cuoco et al. '06ft . This procedure de- 
fines the model, while a single random realization is ob- 
tained by choosing the observed number of events from 
the sources according to the source weights and the 
declination-dependent Auger experimental exposure. 

The model thus depends directly only on n s and the 
choice between sources distributed uniformly or accord- 
ing to the PSCz catalogue (of course, implicitly it also 
depends on the assumptions of sufficiently small mag- 
netic field deflections). Additionally, the model depends 
via the weights of Eq. ([2]) on the type of primary parti- 
cle, the energy spectrum of the sources, and the energy 
scale and cut E cut . The latter dependence arises, be- 
cause the energy scale of CR experiments has a relatively 
large systematic error that is difficult to d etermine. In 
particular, it has been argued ([Teshima '071 ) that the en- 
ergy scale of Auger should be shifted up by 30-40% to 
obtain agreement wit h the spectral shape predicted by 
e + e~ pair production ([Berezinskv et al. '061 ) and the CR 
flux measured by other experiments. Therefore wc use 
two different values for the energy cut, E cut = 60EcV 
assuming that the energy calibration of Auger is cor- 
rect and -Ecut = 80EcV inspired by the dip interpre- 
tation. We do not include in this work the finite en- 
ergy resolution of the Auger experiment that is of or- 
der 20% in AE/E. A finite energy resolution would re- 
sult in an effective decrease of the nominal energy cut 
due to the ste eply falling CR spectr um and to a larger 
GZK horizon (Kachelr iefi et al. '07ft . Similarly, we do 
not account for the stochasticity of the energy loss in 
the photo-pion production process. Both of the effects 
are subdominant at the moment with the present low 
statistics while a more careful analysis will be needed 
in the future when more data is available. For a rough 
estimate of the influence of both effects on n s one may 
compare how the 30% up-ward shift of E cn t from 60 to 
80EcV changes n s . Finally, throughout this work we 
consider proton primaries, but note that the combina- 
tion of nuclei with l arge deflections and few sources has 
been advocated too (jArmengaud et al. '051 : iFargion 'Oil : 
IGorbunov et al. '08[ ). A further brief discussion of this 



point is postponed to Sec. 14.21 

The cumulative autocorrelation of the data C*($) , 
which is a single, one-dimensional function, has now to be 
compared with the hypothesis X , for which we have vari- 
ous Monte Carlo realizations Ck{&\X), k = 1, . . . , M, (we 
use typically M ~ 10 5 ). A standard method to compare 
data and model is to use angular bins §i so that to substi- 
tute the continuous function C(#) with the discrete set of 
values C(#j). The Monte Carlo realizations can then be 
used to calculate the marginalized probability distribu- 
tion of each single C(&i\X) or, if required, the joint prob- 
ability distribution of two C variables C('di\X),C('&j\X), 
three C variables or more. In practice, to derive the best 
fit value n s as well as the goodness-of-fit for the cho- 
sen hypothesis X a possible way is to calculate the mean 
(C($i|X)} and the variance Oi per bin as well as the corre- 
lation matrix Uij and then to perform a % 2 test. However, 
the difficulty to deal with such an high dimensional prob- 
ability space and the generally strong non-gaussianity of 
the probability distribution make the x 2 method clearly 
not optimal for the problem at hand. 

The usual way to circumvent the problem is to use 
the Monte Carlo to calculate the chance probability to 
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Fig. 1. — Penalized chance probabilities p_(n s ), p+(n s ) and p(n s ), for £? C ut = 60 EeV (top panels) and E C ut = 80 EeV (bottom panels). 
The left column reports the case for uniformly distributed sources, the right panel for sources following LSS with the bias of the PSCz 
Galaxy catalogue. Also shown is the 95% and 99% confidence level. 



observe stronger clustering than in the data. Given the 
problem at hand we slightly generalize the method defin- 
ing two functions: the chance probability to observe 
stronger (P + (d\X)) or weaker (P_(i9|A)) clustering at 
the angular scale i? than in the data given by, 

1 M 

p±w) = jjYs^Mamx) - am . (3) 

i=l 

The minimum of the chance probability can then be used 
as a global estimator for the agreement of the hypothesis 
with the data. In our case with a scan of P±('d\X) over 
$ we obtain the two minima P_(Jf) for the minimum of 
P-(&\X) and P+(X) for the minimum of P + (&\X), re- 
spectively. The simplest way to combine these two esti- 
mators is then to use the product P(X) = P- (X)P + (X). 

However, a drawback of this method is that nei- 
ther the quantity P(X) or the single P±(X) are 
truly probabilities. The scan over the angular 
scales fl, in fact, introduces a bias which need 
to be corrected with the use of a penalty fac- 
tor dTinvakov fc Tkachev 'Oil iTinvakov fc Tkachev '03l : 
iFinlev fc Wester hoff '04f ). More precisely, to obtain cor- 
rect probabilities the identical procedure as described 
above needs to be performed with many simulated data 
sets. Counting how often smaller values of P±(X) and 
P(X) are obtained by chance with respect to the case 



of the data, provides true penalized chance probabilities 
p+PQ, p_(X), andp(X). 

The use of the chance probability tool has often gener- 
ated some confusion in the past. An emblematic case 
is the significance of the small scale clustering in the 
AGASA data for which very different estimates ranging 
from ~ 10~ 6 to ~ 10~ 2 ha s been reported in various stud- 
ies (see for example Rcf. (|Takcda et al. '99l )) depending 
on the use or not of the penalty and on different a priori 
choices of the angular an d energy scale of reference (see 
(jFinlev fc Westerhoff *04l ) for a detailed account). The 
effect of the scan can then be quite relevant and a major 
point in the following is that the effect of the penalty 
is correctly taken into account when quoting the con- 
straints on n s and the constraints are further compared 
with the case of an a priori choice of the relevant angular 
scale. 

Note, anyway, that the penalty calculation can be 
avoided if a particular angular scale is chosen a priori 
and the values of the chance probability at this scale 
are employed. However, the scan over all angles avoids 
possible bias, in contrast to the choice of a single an- 
gular scale, which introduces some theoretical prejudice 
even if this choice may be physically motivated. In the 
case at hand, it is unclear if this should be dominated 
by: The ~ 1° angula r resolution of the detector, as 
in (jTakami fc S ato '08) which implicitly assumes neg- 
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ligiblc magnetic field deflections; by a ~ 3° scale, as 
suggested by the cross-correlatio n with active galact ic 
nuclei (AGN) revealed by Auger (|Abraham et al. 'OTP 1 ; 
or yet some other scale, as the 6° separation considered 
in ( Abraham et al. '08d ). To emphasize this dependence, 
we summarize the results of this kind of "single-bin" 
analyses in table [T] and compare them with the ones ob- 
tained with the global method, reported in the last col- 
umn (see next Section). Note, in p articular, how the 6° 
bin chosen in ((Abraham et al. '08cD leads systematically 
to an overly stringent bound. 

3. INTERPRETATION 

In Fig.QJ we report the results for the quantities pi(X) 
defined above, where X = {n s , E cut , k}, the latter two 
being in our case two-valued discrete variables. Because 
the number of CRs usable for this analysis is still small, 
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Fig. 2. — Cumulative autocorrelation function for iYAueer =27 
events compared with the best fit model for the uniform (top) and 
LSS (bottom) cases. In the latter case, only 22 events survive 
outside the PSCz mask. Both plots assume with E cu t = 80 EeV. 
The mean model autocorrelation is shown together with the la 
and 3<r regions. 

1 Note that, for un-correlated deflections, the window size to use 
for autocorrelation studies would be v2~ X 3.1° ~ 4.3°; actually, 
since deflections from the source are correlated and the energies of 
events similar, the relative deflections for a single chemical species 
would be likely < 3°. 



all four hypotheses are compatible with the data at the 
2cr level for some range of n s values. Yet, several inter- 
esting conclusions can already be drawn. The best fit is 
achieved for sources following the PSCz distribution and 
a source density n s ~ 1 x 10 _4 /Mpc 3 . Also, indepen- 
dently of -E cu t, both the penalized probability and the 
range of n s with compatibility at 95% C.L. arc larger for 
the LSS model than for the uniform case. Reversing the 
argument, as can be read more quantitatively in table [1] 
we can see that the constraints are generally stronger for 
the uniform cases with respect to the LSS ones, but this 
is achieved only at the expense of a worse general best 
fit. The fact that the LSS models fit better the data is 
not surprising: Most of the Auger event are aligned along 
the local overdensity known as the Supcrgalactic plane 
which is suitably reproduced with the use of the PSCz 
catalogue within our LSS scenario. 

The case of a uniform distribution of "infinitely many" 
sources (n s — > oo each with an infinitesimal luminos- 
ity), is excluded for both energy cuts at the 95% C.L.: 
The upper bound is n s < (1 3) x 10~ 4 Mpc -3 . This 
is another way to say that the Auger data are inconsis- 
tent with a structureless UHECR sky, independently of 
the use of a catalogue and of a pre- determined angular 
scale for the search. This is, in our opinion, an impor- 
tant milestone in the development of UHECR astron- 
omy. While the best fit point for n s is approximately 
a factor 10 higher than found in earlier studies using 
AGASA da ta above E cnt > 40 EeV (in the AGASA en- 
ergy scale) (lYoshigushi et al. 'OllBlasi fc De Marco '041 ; 
iKachelriefi fc Semikoz '05f h the shape of the chance 
probability p(n s ) agrees: For low values of n s , p(n s ) is a 
steeply decreasing function of n s , since the probability to 
observe multiplets from the same source increases fast. 
In particular, the radius within which 70% of all observed 
UHECRs with energy above E cut = 80 EeV are produced 
is R w 60 Mpc (|Cuoco et al. '06D 2 . As a result, the num- 
ber of sources within this radius becomes less than the 
number of observed CRs events for densities smaller than 
n s « 10~ 5 . Such a scenario would require large deflec- 
tions (and probably nuclei primary) and thus contradicts 
our assumptions. On the other hand, p(n s ) decreases 
relatively slowly for high densities and only weak con- 
straints can be obtained with the current data set for the 
maximally allowed value of n s . Since both an increase of 
E cut and of the bias of the sources leads to a decrease of 
the effective number of sources inside the GZK volume, 
large values of n s have the strongest constraint in the 
case of uniformly distributed sources and E cut = 60 EeV 
(left, top panel) and weakest for sources following the 
LSS and E cnt = 80 EeV (right, bottom panel). 

In Fig. [5] we show the model autocorrelation function 
with la and 3c shaded regions for the best fit uniform 
and LSS model for E cut = 80 EeV, both corresponding 
to n s « 1.4 x 10 _4 Mpc~ 3 , together with the data. At 
small angular scales, ■& < 3°, the data show a deficit of 
clusters compared to the expectation for the the best 

2 The quoted value depends on the use of the continuous- 
energy loss approximation, the actual value increasing to ss 70 Mpc 
due to the stochastic nature of the photo-pion production 
and to 100 Mpc furth er considering a 20% energy resolu- 
tion IjKachchicfi et al. '071) . For the estimate of n s , however, we 
do not use the concept of horizon size explicitly, which is here in- 
troduced only for illustration. 
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TABLE 1 

The estimated number density of sources (at 95% confidence level) under different assumptions, using only the first bin 
information with different sizes, and compared with the global method. 



fit density from the global analysis. This "tension" is 
qualitatively present in most models fitting the data. 
Within our assumptions, this deficit is explained in a 
natural way by (relative) deflections of this size in mag- 
netic fie lds. This value is com parable with the 3.2° found 
in Ref. ([Abraham et al. '08cD that optimizes the correla- 
tions of the same data set with AGNs. The absence of 
small-scale clusters is also responsible for the shift in the 
best fit value for n s compared to old analyses using the 
AG AS A data. The result thus shows that, intriguingly, 
also the autocorrelation function can be employed as a 
sensitive tool for magnetic field studies. Clearly, how- 
ever, a more detailed study needs to be complemented 
with a model of the intervening magnetic fields while, 
likely, more statistics is required to derive significant con- 
straints. 

Our result is in contrast with the one of 
(jTakami fc Sato 'Oct ) which instead report a small 
scale clustering in the Auger data. Notice how- 
ever that, with res pect to our work, the authors of 
([Takami fe Sato '611 ). although including an explicit 
treatment of galactic and extragalactic magnetic fields, 
make use of the non-cumulative autocorrelation func- 
tions and do not take into account penalties for the 
scan over the angular scale. Their claim of small scale 
clustering within 1° is thus likely to disappear if the 
angular scan is taken into account and the comparison 
is made properly with respect to the a mod el effectivel y 
fitting the data. Thus, as already noted in (jHarari '04h . 
the interpretation of small scale clustering could be 
misleading if it is not defined with respect to a proper 
model of the distribution of sources. Apart from this 
point, however, the constr aints we obtain on n s are 
roughly in agreement with (jTakami fc Sato '"08f ). This 
is mainly due to the very low statistics available at 
present which is not still sensitive to the exact analysis 
method employed . We stress however, as noted in 
(|Cuoco et al. '08bl ). that already with a statistic of ~3 
times the present one a formally correct analysis will 
become crucial to avoid biased results. 

Figure [2] also clarifies why the LSS case gives a better 
fit to the data: As can be seen, the LSS best fit model fits 
nicely the data, basically within ler over all the angular 
scales, while the best fit case for uniform sources shows a 
~ 2cr deficit in the broad range 4°-30°. A lower n s (i.e. 
an higher clustering) would not help because it would 
give much more pairs than the data in the l°-4° range. 
Thus, at the end, even with the best possible compromise 
the agreement with the data is only at the 20% level 
for the best fit (or, equivalently, the uniform model is 
excluded at the 80% C.L.). 

We summarize in table[T]the list of best fit n s with 95% 
error bars for the four cases considered and for different 
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Fig. 3. — Contours of equal chance probability p(0|n s ) for the 
LSS case and E C ut = 80 EeV. 

choices of the angular scale or for the global autocorre- 
lation analysis. The crucial point to notice here is that 
the derived n s intervals are sensibly biased with respect 
to each other for different choices of the a priori angu- 
lar scale. Thus, the choice of the angular scale crucially 
affects the result and should be avoided unless the given 
scale has some strong physical motivation. This prob- 
lem is of course avoided employing a global comparison. 
The choice of a single angular scale also affects the error 
bars which can be both larger or smaller with respect to 
the global case. This can easily understood from Fig. [3] 
where we can see that the choice 4°-15° is optimal be- 
cause excess of clustering is observed for low n s while a 
deficit is present for high n s giving thus the tightest con- 
straints. Again, however, a choice in this range is not a 
priori motivated so that from the global analysis we get 
somewhat larger error bars properly taking into account 
all the angular scales. 

Although it is difficult to draw strong conclusions on 
the nature of the UHECR sources at this moment, it 
is worth noting that X-ray selected AGNs with X-ray 
luminosity L > 10 43 erg/s naturally fall in the range 
tiagn = (1 -r 5) X 10- 5 Mpc~ 3 (jSteffen et al. '031 ). at 
the same time, for AGNs with densities in the range 
(10~ 5 -j- 10~ 4 )Mpc -3 , the required luminosity is of 
the order £/n AG N ~ (10 40 + 10 41 )erg/s in UHECRs 
above ~ 60 EeV. So, simple energetics arguments 
are consistent with the inferred values for the den- 
sity; on the other hand, ac celeration efficiency (see 
e.g. (|Ptitsvna fc Troitskv '08f )) tends to favor some 
subclasses of objects and thus a lower inferred density 
of sources, yet typically within the 95% C.L. range 
inferred for n s . Another effect which might reconcile 
the small tension is that the sources can be burst- 
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ing/transicnt or beamed. In this case the luminosity 
requirement is softened (see e.g. the model proposed by 
iFarrar fc Gruzinov '08f ). However, the effectively visible 
number density of sources decreases, too, unless some 
isotropization t akes plac e after acceleration but close to 
the source (see iSigl '081 for more details). Of course, 
these are very preliminary considerations based on a 
limited number of data. For example, these arguments 
might be significantly modified when accounting for 
a more realistic luminosity function. The effects of 
extra-galactic magnetic fields and the consequent UHE- 
CRs time delays can also be quite relevant, although 
at present our knowledge of EGMFs is affected by 
large uncerta i nties and needs t o be better understood 
(see (ISigl '081: iDolag et al. '04l ; iKotera fc Lemoine '08t 
IRvu et al. '08f )). Also, in presence of a heavy-nuclei 
component and/or of extragalactic magnetic fields, a 
significant fraction of events mig ht be associated wit h 
the nearest AGN Cen A (see e.g. iGorbunov et al. '08l ). 
Interestingly signatures in the gamma and neutrino 
bands are expecte d to help di sentangling among many 
scenarios (see (ISigl '08l ; ICuoco fc Hannestad '07 ; 
KachelrieB et al 'OS iBecker fc Biermann 'OS : 

Halzen fc O'Murchadha" 7 ?)! lHardcastle et al. '081 )). 
enlarging the realm of multi-messenger astronomy. 

4. DISCUSSION ON SOME SIMPLIFYING ASSUMPTIONS 
4.1. The role of the energy cut 

The energy-cut used by the PAO to produce 
the sample used in the present analysis was cho- 
sen in order to maximiz e the cross-correlatio n signal, 
see (j Abraham et al. '08bl : lAbraham et al. '08d ) for de- 
tails. One may wonder how this selection affects the 
conclusions of our present work. In principle, the op- 
timal energy cut for the autocorrelation signal and the 
cross-correlation signal with a given catalogue differ one 
from another, although in the case of a common physical 
origin one does expect some correlation between them. 
In particular, it seems reasonable that the optimal cut for 
a global autocorrelation might reside at a lower energy, 
since the small-scale displacement from putative sources 
is not as relevant to the signal as for cross-correlations, 
and the larger statistics helps. Lacking a direct access 
to the data, it is hard to estimate quantitatively how 
large a bias is introduced by focusing on the sample 
of public availabl e events. From the Fig. 2 presented 
in ()Mollerach '07f ). however, one can draw two qualita- 
tive conclusions: i) that a correlation between the two 
optimal cuts in energy is indeed present; ii) that a slightly 
different cut, in the range 40 < E cut /F,eV < 60, should 
still lead to an appreciable clustering of the events. 

We checked that this is indeed the case by performing 
the same analysis as before, but now adding a further 
scan over the range of accessible energy-cuts, i.e. any 
value E > 57EcV. We plot the results in Fig. gp. One 
can note that in all cases the best fit improves and the 
constraints on n s worsen, as expected given the further 
penalty due to the energy scan. Yet, most qualitative 
features described in the previous section stay the same: 

3 Some ripples visible in Fig.[J]are due to the relatively low num- 
ber of Montecarlo: the scan to account for the further penalty fac- 
tor is computationally quite expensive and given the partial nature 
of the answer there is no motivation to refine the results further. 



for example, the LSS model is still preferred over a uni- 
form one. It should be also said that if the true min- 
imum of the chance probability is below E = 57 EeV 
and thus not included in the scan, then these constraints 
are "over-penalized" and thus looser than necessary. At 
the moment, it is impossible to draw more quantitative 
conclusions, since our scan suggests that it is likely that 
the optimal cut for the autocorrelation function is be- 
low E = 57 EeV, a range for which the events are not 
publicly available. 

4.2. Assumptions on the chemical composition 

In this article, we assumed dominant proton primary as 
a basic working hypothesis. It is worth noting, however, 
that the experimental situation on the chemical compo- 
sition at UHE is far from settled: while anisotropy data 
point to a relatively light composition, the results of the 
fluorescence detector of the Pierre Auger C ollaboratio n 
favor a significant fraction of heavy nuclei ()Unger '07f ). 
Yet, for the interpretation of these results one must rely 
on simulations employing hadronic interaction models. 
These are not based on a first-principle theory, rather 
on models calibrated on "low-energy" collider data, then 
extrapolated about two orders of magnitude beyond the 
center of mass energies experimentally probed. 

A proper quantitative assessment on how our conclu- 
sions vary in a mixed composition scenario goes beyond 
the purpose of the present paper. Yet, at a qualitative 
level, we can note that several effects would come into 
play. First of all, in the unrealistic case where one could 
forget about magnetic deflections, the major effect would 
be a reduction of the energy-loss horizon (but for iron, 
whose horizon is similar to the proton one) . This should 
enhance the anisotropy pattern, due to the prominence 
of nearby accelerators. 

When including (the poorly known) magnetic fields, 
two additional effects are relevant: i) for a given en- 
ergy, the higher the charge the larger the deflection and 
hence the loss of information at small angular separa- 
tions. Quantifying how large is this scale is a difficult 
task. We note that protons of these energies in the sole 
Galactic field likely suffer a few degrees absolute deflec- 
tions, see (iKachelriess et al. '05f ) . implying a degree-scale 
smoothing in the relative deflections important for the 
2pcf. This is comparable with the angular resolution of 
the PAO: in this optimal case the whole information in 
the 2pcf starting at i9 m ; n ~ 1° could be used. However, 
a different Galactic field model (especially towards the 
Galactic Center), the presence of heavy nuclei, and/or 
significant extragalactic magnetic fields can easily lift 
by one order of magnitude or more, ii) the other 
effect is that the real path-length of the nucleus would 
be longer than the distance to the source: thus, the dis- 
tance of accessible sources would be even shorter than 
estimated from energy-loss considerations. This is par- 
ticularly relevant if the nucleus spends a lot of time in 
a magnetized region surrounding the accelerator (e.g. a 
magnetized cluster in which it is immersed) before escap- 
ing in the Intergalactic Medium. 

Finally, if the maximal energy of acceleration of dif- 
ferent species of nuclei fall by unfortunate coincidence in 
the same region expected for the GZK feature, slightly 
different energy cuts in the data (as well as statistical and 
systematic errors on the energy scale) might significantly 
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Fig. 4. — Penalized chance probabilities p_(ra s ), p+(n s ) and p(n s ), for i? C ut = 60EeV (top panels) and E C ut = 80EeV (bottom panels) 
taking into account the effect of the energy scan. The left column reports the case for uniformly distributed sources, the right panel for 
sources following LSS with the bias of the PSCz Galaxy catalogue. Also shown is the 95% and 99% confidence level. 



change the expected pattern of anisotropics. The same 
happens if the proportions of different nuclei accelerated 
at the source change as a function of the energy. This 
is in principle a possibility, especially if different classes 
of objects contribute to the events at slightly different 
energies. 

The general pessimistic conclusion is that if several of 
the above effects are relevant (or perhaps a single one 
is dominant) the capability of performing UHECR as- 
tronomy would be greatly reduced. While one still ex- 
pects indications for anisotropics, inverting the problem 
and inferring the source/propagation medium properties 
would require a much larger statistics (especially in the 
trans-GZK region): disentangling the different effects is 
in fact a formidable task. 

To provide a glimpse of how some of the above effects 
alter the reconstruction of n s (our main topic here), in 
Table [2] we report how the constraint on n s degrades 
as a function of a "smoothing angle" i? m i n , below which 
we assume that the 2pcf information is completely lost. 
The main trend is that, if the smallest angular scales 
are neglected, it is easier to find parameter configura- 
tion fitting the data and, correspondingly, the allowed 
range for n s widens. This had to be expected, given the 
shape of the correlation functions shown in Figs. 2-3. 
In particular we find that with t? m i n = 3° we obtain al- 
most the same results as in the global case. The case 



$min = 10° still places useful constraints especially for 
the i? cu t = 80EeV case, while finally using only the in- 
formation above ^ m ; n = 30° basically no constraints on 
n s are obtained. Note however that relative deflection 
angles of that sort would imply overall deflections even 
larger, seriously questioning the perspectives of present 
instruments to perform some form of UHECR astronomy. 

5. SUMMARY 

We used the first UHE data released by Auger to per- 
form a global autocorrelation study of UHECR arrival 
directions, assuming proton primaries. The major ad- 
vantage of our tool is that no biases are introduced by 
the a priori choice of a single angular scale. The main 
observable we have focused on is the number density n s 
of ultrahigh energy cosmic ray (UHECR) sources. While 
the global analysis does not bias n s by what is the the- 
oretical prior of the "relevant" angular scale, still it is 
important to establish how the extraction of the allowed 
range of n s from the data depends on a number of other 
effects, an issue often overlooked in the literature. In 
particular, here we discussed the systematic energy scale 
uncertainty and of the bias of UHECR sources with re- 
spect to Large Scale Structures. As a first attempt to ex- 
tract some information from the data, we compared four 
hypotheses: a structured universe (following the PSCz 
catalogue) and an isotropic case, each for two possible 
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TABLE 2 

The estimated number density of sources (at 95% confidence level) under different assumptions on the minimum angle 
above which the 2pcf information is preserved, 1? m i n . n.c. stands for "no constraints". 



values for the absolute energy scale of the PAO experi- 
ment. The density n s is the free continuous parameter 
in terms of which constraints have been analyzed. 

Not surprisingly, we find that the number of CRs us- 
able for this analysis is still small and all four hypotheses 
arc compatible with the data at the 2a level for some 
range of n s . Yet, several interesting observations can be 
tentatively drawn. The best fit is achieved for sources fol- 
lowing the matter tracer distribution, and a source den- 
sity n s ~ 1 x 10 _4 /Mpc 3 . It is interesting to note that the 
data show some preference (although not significant, yet) 
for a str uctured universe: other recent statistical stud - 
ies, like (jKashti fc Waxman '08l : lKoers fc Tinvakov '09T ). 
qualitatively agree in that respect. Also, there is indica- 
tion that the case of a uniform distribution of "infinitely 
many" sources (n s — > oo each with an infinitesimal lu- 
minosity), is excluded for both energy cuts at the 95% 
C.L.: The upper bound is n s < (1 3) X lO^Mpc -3 . 
This is another way to say that early Auger data sug- 
gest that data are poorly consistent with a structureless 
UHECR sky, independently of the use of a catalogue and 
of a pre-determined angular scale for the search. 

Compared to a benchmark number density of proton 
sources n s ~ 1 x 10~ 4 /Mpc 3 , a factor ~ 2 lower den- 
sities are preferred if the current Auger energy scale is 
correct, a factor ~ 2 higher value if it is underestimated 
as required by the dip model. Including the finite energy 
resolution of the Auger experiment into the analysis will 
reduce further the best fit value for n s . The width of the 
allowed region is dominated at present by the statisti- 
cal error due to the small number of events. Nominally, 
approximately a three times larger sample is needed to 
reduce the Poisson error below the typical differences be- 
tween source candidates. Once that level of statistics is 
reached, other effects will provide the main source of er- 
ror, a major one being the systematics on the energy scale 
as we illustrated here. In the future, other effects such as 
the systematic and statistical errors in the energy deter- 
mination of UHECR events and the stochastic nature of 
the photo-pion process need to be included to correctly 
determine the best fit value of n s . Even considering the 
above limitations, preliminary conclusions arc the follow- 
ing: First, the fit generally improves for sources follow- 
ing the LSS compared to uniformly distributed sources; 
qualitatively, for AGNs which arc known to be even more 
overdensity-biased than the LSS (see e.g. discussion in 
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(|Cuoco et al. '08 rj )). the agreement is expected to im- 
prove even further. In particular, the case of an uniform 
distribution of "infinitely many" sources, i.e. a structure- 
less UHECR sky, is excluded at 99% C.L. Second, the ab- 
sence of clustering on scales smaller than a few degrees is 
most easily understood as the effect of magnetic smear- 
ing of comparable size. This, intriguingly, suggests that 
autocorrelation studies can be employed as a complemen- 
tary tool to study galactic (and extragalactic) magnetic 
fields. Conclusions about the source scenario are compli- 
cated by the the limited knowledge of the EGMFs and 
the presence of beaming and/or bursting effects which 
arc difficult to disentangle with the use of the autocor- 
relation alone. Hence, the true source density could be 
somewhat lower that the best fit value found: The energy 
scale error alone shifts the best fit value by a factor two 
or three; luminosity function effects are likely relevant, 
too. The scarce statistics at the moment is a serious 
limiting factor to constrain more realistic models, but 
with a greater exposure of currently existing instruments 
and the multi-messenger combination of gamma-ray and 
neutrino data, some of these issues will probably be ad- 
dressed and solved in the near future. 

For a significant contamination of nuclei in the sam- 
ple, a much more complicated analysis is needed, since 
many more variables enter the game. Qualitatively, one 
can expect that although future data might allow to dis- 
entangle a proton-dominated sample from a more com- 
plicated picture, inferring source properties and disen- 
tangling them from magnetic effects (in a few words, 
performing UHECR astronomy) should wait for a ma- 
jor jump in exposure, perhaps beyond the capabilities 
of currently planned instruments. Similar considerations 
apply unfortunately to the constraints to the cross sec- 
tions of UHECRs as well. 



ACKNOWLEDGMENTS 

We are grateful to Giintcr Sigl for useful discussions. 
M.K. thanks the Max-Planck-Institut fur Physik in Mu- 
nich for hospitality and support. P.S. is supported by the 
US Department of Energy and by NASA grant NAG5- 
10842. Fermilab is operated by Fermi Research Alliance, 
LLC under Contract No. DE-AC02-07CH11359 with the 
United States Department of Energy. 



J. Abraham et al. [The Pierre Auger Collaboration], Phys. Rev. 

Lett. 100, 211101 (2008b). 
J. Abraham et al. [Pierre Auger Collaboration], Astropart. Phys. 

29, 188 (2008c). 

J. Abraham et al. [Pierre Auger Collaboration], Phys. Rev. Lett. 
101, 061101 (2008d). 



9 



M. Aglietta et al. [Pierre Auger Collaboration], Astropart. Phys. 
27, 244 (2007). 

R. Aloisio and F. Tortorici, Astropart. Phys. 29, 307 (2008). 

E. Armcngaud, G. Sigl and F. Miniati, Phys. Rev. D 72, 043009 
(2005). 

J. K. Becker and P. L. Biermann, arXiv:0805.1498 [astro-ph]. 

V. Berezinsky, A. Z. Gazizov and S. I. Grigorieva, Phys. Rev. D 

74, 043005 (2006). 
P. Blasi and D. De Marco, Astropart. Phys. 20, 559 (2004). 
A. Cuoco et al., JCAP 0601, 009 (2006). 

A. Cuoco, G. Miele and P. D. Serpico, Phys. Lett. B 660, 307 
(2008a). 

A. Cuoco et al., Astrophys. J. 676, 807 (2008b). 

A. Cuoco and S. Hannestad, Phys. Rev. D 78, 023007 (2008). 

K. Dolag, D. Grasso, V. Springel and I. Tkachev, JETP Lett. 79, 

583 (2004) [Pisma Zh. Eksp. Teor. F iz. 79, 719 (2004)]. JCAP 

0501, 009 (2005) astro-ph/0410419] . 

D. Fargion, arXiv:0801.0227] |a stro^pnT. 

G. R. Farrar and A. Gruzinov, larXiv:0802. 10741 [astro-ph] . 

C. B. Finley and S. Westerhoff, Astropart. Phys. 21, 359 (2004). 

D . S. Gorbunov, P. G. Tinyakov, I. I. Tkachev and S. V. Troitsky, 
larXiv:0804.1088l [astro-ph] . 

K. Greisen, Phys. Rev. Lett. 16, 748 (1966). 

F. Halzen and A. O'Murchadha, arXiv:0802. 0887] [astro-ph]. 

D. Harari, S. Mollerach and E. Roulet, JCAP 0405, 010 (2004). 
M. J. Hardcastle, C. C. Cheung, I. J. Feain and L. Stawarz, 

larXiv:0808.1593l [astro-ph] . 
M. Kachelriess, P . D. Serpico and M. Tes hima, Astropart. Phys. 

26, 378 (2006) [arXiv:astro-ph/0510444] . 
M. Kachelriefi and D. Semikoz, Astropart. Phys. 23, 486 (2005). 
M. Kachelriefi and D. Semikoz, Astropart. Phys. 26, 10 (2006). 
M. Kachelriefi, E. Parizot and D. V. Semikoz, larXiv:0711.3635l 

[astro-ph] . 

M. Kachelriefi, S. Ostapchenko and R. Tomas, larXiv:0805.2608l 

[astro-ph] . 

T Kashti and E. Waxman, JCAP 0805, 006 (2008) 
[arXiv:080l!45T6l [astro-ph]] . 



H B. J. Koers a nd P. Tinyakov, JCAP 0904, 003 (2009) 

|arXiv:0812.0860l [astro-ph] ] . 
K Kotera and M Lemoine, Phys. Rev. D 77 (2008) 123003 

|arXiv:0801. 14501 [astro-ph]] . 
S. Mollerach [The Pierre Auger Collaboration], proceedings of the 

"30th International Cosmic Ray Conference", Mcrida, Mexico, 

2007, arXiv:0706.1749] [ astro-ph]. 

K. Ptitsyna & S. Troitsky. larXivT 0808.0367 

D . Ryu, H. Kan g, J. Cho and S. Das, Science 320, 909, 

arXiv:0805.2466 [astro-ph]. 
W. Sau nders et al, Mon Not. Roy. Astron. Soc. 317, 55 (2000). 
G. Sigl, arXiv:0803.3800 [astro-ph]. 

G. Sigl, F. Miniati and T. Enfilin, Phys. Rev. D 68, 043002 (2003) 
Phys. Rev. D 70, 043007 (2004). 

A. T. Steffen et al. Astr ophys. J. 596, L2 3 (2003). 

H. Takami and K. Sato, [arXiv:0807.3442 [astr o-ph] . 

M. Takeda et al, Astophys. J. 522, 225 (1999) |astro-p h/9902239 . 

See also: P. G. Tinyakov and I. I. Tkachev, JK1F Lett. 

74, 1 (2001) [Pi sma Zh. Eksp. Teor. Fiz. 74, 3 (2001)] 

astro-ph/010210ll; D. De Marco, P. Blasi and A. V. Olinto, 

JCAP 0601 (2006) 002 [arXiv:astr<>p h/0507324 ; H. Yoshiguchi, 

S. Nagataki and K. Sato, Astrophys. J. 614, 43 (2004) 

[astro-ph/0404411] . 
M. Teshima, rapporteur talk at the "30th International Cosmic 

Ray Conference" , Merida, Mexico, 2007. 
P. G. Tinyakov and 1. 1. Tkachev, JETP Lett. 74, 445 (2001) [Pism a 

Zh. Eksp. Teor. Fiz. 74, 499 (2001)] [arXiv:astro-ph/0102476] . 
P. Tinyakov and I. Tka chev, Phys. Rev. D 69, 128301 (2004) 

arXiv:astro-ph/0301336 . 

M. Unger [The Fierre Auger Collaboration], larXiv:0706.1495l 

[astro-ph] . 

H. Yoshiguchi, S. Nagataki, S. Tsubaki and K. Sato, Astrophys. J. 

586, 1211 (2003) [Erratum-ibid. 601, 592 (2004)]. 
G. T. Zatsepin and V. A. Kuzmin, JETP Lett. 4, 78 (1966) [Pisma 

Zh. Eksp. Teor. Fiz. 4, 114 (1966)]. 



